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An implementation of premixed equilibrium chemistry has been completed 
for the OVERFLOW code, a chimera capable, complex geometry flow code widely 
used to predict transonic flowfields. The implementation builds on the computa- 
tional efficiency and geometric generality of the solver. 


x, y , z cartesian position coordinates 


Introduction 

Building on the successful application of the 
perfect gas version of OVERFLOW to hypersonic 
flowfields, 1 the next step in implementing finite 
rate chemistry capability in OVERFLOW is the 
generalization of the gas model from single species 
perfect gas to allow pressure and temperature as 
arbitrary functions of density and internal energy. 
The specific thermodynamic model used is that of 
Liu and Vinokur. 2 This model provides pressure 
and temperature for an equilibrium mixture of air 
over a wide temperature and pressure range. 
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Method 

The Reynolds averaged Navier-Stokes equa- 
tions, written in conservation law form are 


9Q 0F 0G 0H: 
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Nomenclature 

y pc 2 /p 

e internal energy per unit mass (/Jc v dT) 

£r total internal energy per unit mass e + 
j(u 2 + v 2 + w 2 ) 

e internal energy per unit volume (J Q pc v dT) 

k turbulent kinetic energy j Y u'u' 

Pt turbulent eddy viscosity 

p mass per unit volume 

0 angle along sphere surface (degrees) 

c isentropic sound speed 

hj stagnation enthalpy (e + p/p 4- \ (u 2 +V 2 + 

w 2 )) 

L reference length 

p static pressure 

Q conservative variables [p, pu, pv, pw, pet] 
s entropy 

Sy Strain rate tensor j j 

T static temperature 

u, v, w cartesian velocity components 
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Where 
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[p, pu, pv, pw, pe] (2) 

pu, pu 2 + p + t xx , puv + T X y , puw + r xz , 
pu( e T + P + Ax) + vT xy + WT XZ + q x ] (3) 

PV, PVU + Ty X , PV 2 + P + Tyy , PVW + Ty Z , 
pv(e T +P+Tyy)+WT yz + UT:y x -|-qy] (4) 


|I| we need to be able to get the sound speed 
in terms of these variables and their derivatives. 

From the Gibbs equation for a simple compress- 
ible substance, 


[pu, , pwu + T zx , pwv + Tzy , pw +v + T zz , 
pw(e T + P + t zz ) -I- UT ZX + vT Z y + q z ] (5) 


where, we assume a simple Boussinesq type stress- 
strain relationship for the turbulence, 

Ty = 2(|x 4- M-t ) ( S i j — ^S nn ) — -pKSq (6) 

and the heat flux is given by the Fourier relation 

„ i M- , Ft , 3h 
qi = + <7) 

The thermodynamic relation for pressure p = 
p(p, e) is used in place of the perfect gas assump- 
tion (p = p(y— 1 )e) in computing the fluxes of the 
momentum and energy equations, and the ther- 
modynamic relation for temperature T = T(p,e) 
enters into these equations indirectly, in the eval- 
uation of molecular viscosity and thermal conduc- 
tivity. 

Implicit solution requires the calculation of Eu- 
ler Jacobians (f^-, , and ^=f) and for diag- 

onalized methods, their eigenvectors and eigen- 
values. The Euler fluxes are simply the Navier- 
Stokes fluxes with p., p.j> and k, set to zero. This 
generally involves an evaluation of the isentropic 
sound speed and derivatives of pressure with re- 
spect to p and pe = e. The sound speed is derived 
from the partial derivatives of the pressure re- 
turned by the thermodynamics module, and its 
current implementation allows a capability more 
general than mixtures of perfect gases in chemical 
equilibrium. 

This solver is actually capable of predicting 
flows of any simple compressible substance where 
pressure and temperature are functions of only 
density and internal energy. The sound speed is 
given by the general relation for isentropic sound 
speed. c 2 = As we have a table lookup rou- 
tine which returns values of p, T, |^, f| 


ds = 


we get the relations 


(de +pdv) 


de — pdp/p" 


= P/P 


= 1/T 


which, along with the trivial relations 


allows us to transform the partial derivatives we 
have naturally as a result of the bicubic spline: 
, I 2 into derivatives where [p,s] are the 

° p e 0G p 

independent variables. Thus we may calculate the 
sound speed from 


c 2 = 


3p 9e 
3e p 3p 


= S* + ~2r < n > 

9p € P- 0e p 

The derivatives of pressure required in the 
eigevector determination are 


e 0p 
P 3e P 


(note the distinction between e, as opposed to e). 
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Implementation 

The base code used is OVERFLOW 3 * which 
is a finite difference, chimera capable, complex 
geometry flow code -widely used for perfect gas 
prediction. Matrix dissipation is used with cen- 
tral space differencing and multigrid to provide 
a rapid, robust solution method suitable for high 
speed flow with strong shock waves. 1 

OVERFLOW was modified to replace the per- 
fect gas model built into the code with an arbi- 
trary gas model, where the pressure and sound 
speeds are functions of the density and internal 
energy, as given in Liu. 2,4 An efficient table 
lookup method returns the pressure and tempera- 
ture and their derivatives corresponding to a given 
value of [p, e], A bicubic spline is used to produce 
smooth interpolations of both functions and their 
derivatives. 

Four additional variables are added to the con- 
servative variables [p, pu, pv, pw, pet]. The four 
additional field variables axe fp/poo, X, Pr, T/ToJ. 
These variables are obtained from the conserva- 
tive variables at the start of each iteration, before 
the laminar viscosity is calculated. The perfect 
gas assumptions are then replaced by the use of 
the calculated pressures (in the momentum and 
energy fluxes, and their derivatives), tempera- 
tures (in the viscosity and fconductivity calcula- 
tions), and sound speed (in the implicit Jacobian 
calculations). 

2nd-order central flux differencing combined 
with matrix dissipation 5, 6 is used in forming the 
right hand side residuals, similar to the scheme 
used for perfect gas solutions. 1 The implicit solu- 
tion method used is the Pulliam-Chaussee diag- 
onal algorithm, with the eigensystem 7 the same 
as used to implement the matrix dissipation. Roe 
averaging is used in forming the matrices in the 
dissipation scheme. The pressure used to form 
the pressure switch function is that given by the 
p = p(p, e), as opposed to the perfect gas assump- 
tion pressure. 

The dissipation settings are similar to those 
suggested previously, 1 with the choices (unless 



a) Pressure 


b) Temper- 
ature 


c) Density d) h t on Grid 


Fig. 1 Mqo = 32, p = lOPa = Blunt Body Flow- 
field 


otherwise noted) as 


dis 2 

= 4 

dis 4 

= 0.2 

v ei 

= 0.3 

V £n 

= 0.3 


(14) 


The settings for the initial (coarse grid) flow setup 
have modified eigenvalue limiters, (V e , = V 6n = 
1 .0), which allows the shock to move off the body 
surface into the flowfield. 

The multigrid method 6 has been modified to 
enhance the robustness of the method. The orig- 
inal multigrid method used a linear interpolation 
of the conservative variables to map the coarse 
grid solution onto the next finer grid. This inter- 
polation is replaced with a linear interpolation of 
the variables [p,u,v,w, e], which ensures a posi- 
tive value of e on the finer grid if e is positive on 
the coarser grid. 


Results 

Inviscid Blunt Body Flow 

This case is the inviscid flow past a blunt cylin- 
der (Fig. 1) over a range of freestream static 
pressure and Mach number. The pressure con- 
ditions range from p = IPa (h=80km) to p = 
10KPa(h=16km). The Mach numbers range from 
Moo - S (2.4km/sec) to = 32 (9.5km/sec). 
The solutions were obtained on a 101 x 129 grid 
suitable for computation of a perfect gas, so much 
of the grid lies outside the shock. 
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Fig. 3 Inviscid Cylinder Blunt Body Stagnation Line State 

These solutions demonstrate the ability of the fer predictions in viscous flowfields. The pressure 

matrix dissipation to capture the shock over the profiles in all cases are collapsed nicely when nor- 

wide range of conditions representative of earth malized by the Newtonian pressure, pu^,. 

atmospheric reentry. The smoothing parameters 

, , . r , . Laminar Flat Plate 

are unchanged over this entire range of condi- 
tions, and only the M = 32 cases required an The ability to correctly predict heat transfer 
adjustment of the CFL(halved from the other con- and skin friction is as important as inviscid shock 
ditions) to obtain convergence (Fig. 2). capturing capability. To check the viscous predic- 

While all the cases show a similarity in the tion capabilities, solutions for a laminar flat plat 

stagnation pressure profiles (Fig. 3), the vari- at freestream Mach numbers of 2 and 4 were ob- 

ations of the temperature and density are more tained on 3 grids, 129 x 129, 65 x 65, and 33 x 33. 

pronounced as the Mach number increases. At a The finest grid had a wall normal grid spacing 

fixed Mach number, a lower freestream pressure of 2 x 10 6 (in units of plate length), and grid 

implies a lower post-shock temperature, as the stretching factors less than 1.05. The coarser 

dissociation is more complete at lower pressure. grids are obtained from the finer grid by bicima- 

An interesting trend is evident in the density pro- ti° n (deleting every other point). Uniform inflow 

files. As the Mach number increase from 8 to 24, conditions are applied at the upstream and outer 

the post shock density increases, but from Mach boundary. Simple extrapolation is used as the 

24 to Mach 32, there is a decrease in the post boundary condition at the outflow plane, and the 

shock density especially at the higher pressures. constant temperature no slip is specified on the 
The preservation of total enthalpy across the plate itself. A short region (12, 6 and 3 points 

shock is good at all Mach numbers; Along with respectively on the fine, medium and coarse grid) 

being a check on the fidelity of the solution, this specified as inviscid wall upstream of the plate, 

will be important in obtaining good heat trans- The Reynold number of the flow (based on 
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Fig. 4 La m inar Flat Plate Skin Friction Predictions 



Re, Re, 


a) M .00 — 2 b) M-oo = 4 

Fig. 5 Laminar Flat Plate Heat Tr ans fer Predictions 

plate length) is 600K for = 2, and 300K for Turbulent Flat Plate 
Mqo = 4. The skin friction predictions are all 

nicely in agreement with theory (Fig. 4) and the As turbulent flow is another of interest 
heat transfer predictions (Fig. 5) are similarly for users of the code ’ a turbulent flat P lat R e case 
good. These predictions are obtained with the was com P uted ' Tbis is a = 5 case, 8 with 
same dissipation parameters used to produce the “ adiabatic waib A s P ecial version of tbe ^ 
inviscid wedge and cylinder solutions. Also shown P ro P erties table was created for tbis ca5e > “ tbe 

are the skin friction predictions with the lowered waU tem P erature is approximately 100°K lower 

dissipation settings, which agree quite well with wben the actual variation of s P ecific heat is used ' 
the ‘standard’ dissipation settings. - The modlfied Properties table had; constant spe- 

cific heats and Prandtl Number, with y = 1 .4 and 
Pr = 0.73, and Pr t was chosen as 0.88, consistent 
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c) Moo = 24 d) Moo = 32 

Fig. 2 Inviscid Cylinder Blunt Body Residual 
Histories 



Fig. 6 Turbulent Flat Plate Grid (every 4th 
grid line shown) 

with the earlier computations. 8 The Reynolds 
number based on plate length is as 140 Million, 
with a freestream static temperature of 273.15°K. 

The solution is obtained on two grids, a fine 
grid of 129 x 257 (streamsize x wallnormal), and 
65 x 129, the coarser grid as usual created from 
the finer by removing every other point. The 
grid is created by exponential stretching in both 
streamwise and wall normal directions with ini- 
tial spacings of 1 CP 4 and 1 0 -6 of the plate length 
in each direction, respectively for the fine grid. 
The maximum stretching ratio of the the grid 
spacing is 1.05 streamwise and 1.04 wall normal. 
The outer boundary of the grid contains the weak 
shock created by nose of the flat plate. The grid 
wall normal spacing has Ay^ < 1 for. the coarser 
grid, with AyJ < 0.1 at the end of the fine grid. 

A comparison of predictions using, the 





Fig. 7 Turbulent Flat Plate Skin Friction 


Baldwin-Lomax(BL) algebraic model, the 
Spalart-Almaras 9 (SA) one equation model, and 
the two equation Menter Shear Stress Trans- 
port 10 (SST) are shown in Fig. 7. The SA 
implementation used is the published version, 
as opposed to the modified version standard in 
OVERFLOW. The skin friction predictions ob- 
tained are in good agreement with the Van Driest 
Theory, consistent with earlier work. Again, 
as with the laminar flat plate, solutions were 
obtained with the smoothing parameters halved, 
and the skin friction prediction predictions are 
again insensitive to the smoothing parameter 
choice(Fig. 8). 

Viscous Sphere 

This case considers a laminar viscous flow 
over a 1m radius sphere at Moo = 17.6, p = 
57.4Pa, Too = 200° K, as shown in figure 9. 
This case is the benchmark case for LAURA, 
(http://hefss.larc.nasa.gov). The viscous wall 
boundary condition is no slip, T w = 500° K. The 
finest grid used is 65 (streamwise) by 129 (wall 
normal), with a constant wall grid spacing of 
0.5 p.m. Two coarser grids were derived by this 
one by choosing every other point (33 x 65), and 
every fourth point (17 x 33). The solution is also 
obtained on a full 3D grid, with grid dimensions of 
33 x 65 x 129 (Streamwise x Circumferential x Wall 
Normal). The solution for the 3D grid lies on top 
of solution of the fine axisymmetric grid, and is 
both a check on the -axisymmetric option and the 
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Fig. 8 Smoothing Parameters Sensitivity - SA 
Model 


3D capability of the solver. 

The stagnation profiles for this case (Fig. 10) 
are given in both external flow (linear in axial 
distance) and boundary layer(logarithmic in axial 
distance) coordinates. The pressure profile looks 
similiar to the inviscid case (Fig. 3), and agrees 
well with the prediction of LAURA, but here the 
stagnation enthalpy and static temperature varia- 
tions axe dominated by that of the boundary layer 
variations. The temperature predictions show the 
same boundary layer “lumpiness” predicted by 
LAURA, whereas the total enthalpy prediction is 
quite smooth over the same region. The 3D so- 
lution prediction is denoted by a dotted fine in 
this figure, and follows the fine grid axisymmetric 
solution. Heat transfer predictions(Fig. 11) are 
similiar to those of the LAURA code. 



Fig. 9 Viscous = 17.6, p = 57.4Pa 1m 
Sphere 




Fig. 10 Stagnation Streamline Profiles, — 
17.6, p =57.4Pa lm Sphere 
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Fig. 13 Wall Pressure Distribution, T10-917 
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High Enthalpy Nozzle Expansion 

This case highlights the ability to predict flow- 
fields with extreme expansion. T10-97 is one of 
the First Europe-U.S. High Speed Flow Database 
test cases. 11 This case documents the flowfield in 
the ONERA F4 high enthalpy facility. The re- 
sevoir conditions are p = 373bar, ho = MOT^, 
and the exhaust conditions are Too = 273° K, 
Poo — 68Pa. The resevoir is charged with syn- 
thetic air, composed of |C >2 and |N 2 by mass 
fractions. 

A sequence of three grids is again uti- 
lized, with the finest grid 513(Streamwise) x 
1 29( WallNormal) , and the coarser grids again cre- 
ated via ’bicimation’ (removing every other point). 
The inflow plane is given p t = constant, h t = 
constant boundary conditions, the walls are no 
slip, fixed temperature (T w = 300° K), and the 


Fig. 14 X33 Surface Grid System 

outflow plane is 0 order extrapolation. 

The solutions for this case depend critically on 
the air thermodynanmic models, 11 and the cur- 
rent results (Fig. 13) are consistent with the 
results of other equilibrium air codes utilized to 
simulate this case. The predicted exit velocity is 
4.68km/sec, consistent again with predictions by 
other equilibrium codes. 11 

X33 

A final case is the X33, as solved in the pre- 
vious perfect gas work. 1 This is a multiple zone 
case, at a freestream Mach number of 5.98, an- 
gle of attack of 40°, and Reynolds number of 3.3 
Million. ’Wall temperature is constant at -300° K. 
The GASP and' ’perfect gas’ solutions are taken 
from that earlier, paper, as is the grid system and 
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Fig. 15 Surface Pressure, X33 
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Fig. 16 Centerline Pressure Comparison 
boundary conditions. 

The surface pressure axe in good agreement 
with the previous perfect gas and GASP solu- 
tions. 12 Surface heat transfer on the windward 
side of the vehicle is somewhat lower than that of 
the experiment, 13-15 where the perfect gas solu- 
tion is higher. This is consistent with what one 
would expect from finite rate chemical reaction 
effects. 


Discussion 

The conclusions of the previous paper 1 for per- 
fect gas cases continue to hold for the general- 
ization of the thermodynamic model of overflow. 
Matrix dissipation yeilds predictions of similar 
quality as upwind MUSCL schemes, at much 
lower computational cost. Here the case can be - 


1 


0.5 


0 

0 0.5 x/L 1 

Fig. 17 Centerline Heat Transfer Comparison 

made for Mach numbers up to 32, and over al- 
titudes from sea level to 65km. A single set of 
parameters for the smoothing worked cases from 
inviscid cylinders to laminar and turbulent flat 
plates, as well as 3D geometries such as a sphere 
and the X33, and for extreme expansions such as 
the ONERA F-4 nozzle flow. The grids used were 
not “special” in any sense, and were constructed 
using standard tools (or re-using grids from other 
solvers), and the starting method worked reliably 
over all cases tried. 

Multigrid had the same beneficial effect as seen 
in the perfect gas solutions, improving robustness, 
and accellerating convergence. Grid sequencing 
and multigrid work were utilized in every solution 
in this paper, and greatly aid the solution process. 

Another piece of transonic technology which 
continued to work well with equilibrium chem- 
istry is the Pulliam-Chausee diagonal algorithm. 
Although coupling this solution algorithm with 
chemistry is problematic, for solutions where the 
mean flow and chemistry systems can be solved 
in an uncoupled m ann er, this algorithm promises 
large computational cost savings. 

The solution of the X33 system shows the feasi- 
bility of using the chimera system for high speed 
flowfields. The ability to use chimera to model 
complex geometries greatly facilitates the ability 
to obtain predictions for flight vehicles, and al- 
lows the tailoring of grid systems component by 
component. For rapid evaluation of an evolving 
system design, this capability is extremely impor- 
tant. 


E Experiment 
« OVERFLOW(PG) 
o GASP 

• OVERFLOW(EQ) 
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Conclusions 

The present method shows promise for rapid 
prediction of complex flowfields where fully mixed 
chemical equilibrium is an appropriate approxi- 
mation, and is another step in the overall goal 
of adding chemically reacting flow capability to 
OVERFLOW. 
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